#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<time.h>
#include<direct.h>
#include<sys\\stat.h>
#include<sys\\types.h>

#define mSol 1
#define mTerra 3.003167261575587e-6
#define mJupiter 9.543512141169373e-4

#define G 4*3.141592653589793*3.141592653589793
#define perielioTerra 0.983235919814466
#define vPerielioTerra 6.389661159796174
#define excentricidadeTerra 0.0167

void inicializacao(void);
void keplerT(void);
void criarArquivos(void);
void F(double *dx, double *z, double *distancia, double m,double M, int dimensao);
void calcularDistancia(double *z,double *distancia, int dimensao, int corpos);
void jacobianaK(double* j, double M,double m,double rx,double ry,int tamanho);
void produtoEscalarVetor(double* v,double a,int tamanho, double *resultado);
void produtoEscalar(double* v,double* v2,int tamanho, double *resultado);
void somaVetores(double* v,double* v2,int tamanho, double *resultado);
void produtoMatrizVetor(double *J,double* delta,int tamanho, double *resultado);
void escreverArquivos(void);
void imprimir(double* v,int tamanho);
double calculaNorma(double *v,int dimensao);
void criarMatrizRK45(void);
void finalizar(void);

double *K1,*K2,*K3,*K4,*K5,*K6;
double  e,a2,m,M,tFinal,dt,*a,*b4,*b5,*c2,*y,*Y,*norma,*temp,*temp2;
double  *dx, *distancia,tempo,*J,*z,*z1,*z2,tempo;
int dimensao, n, T, corpos,escrever,tamanho,renormalizacao,orbitas=1,stages=6;
char *exemplo,*aux,*SO;
clock_t inicio,fim,total;
FILE *f[9];

int main(int argc, char *argv[])
{
	inicio=clock();

	SO=argv[1];
	escrever=atoi(argv[2]);
	tFinal=atof(argv[3]);
	n=atof(argv[4]);

	inicializacao();

	int i;
	/*for(i=0;i<12;i++)
	{
		z[i]+=0.0001;
	}*/

	criarArquivos();
	escreverArquivos();
	criarMatrizRK45();

	//RK4
	int j;
	for(i=1;i<T;i++)
	{
		//temp2=soma parcial

		//RK4
		//stage 1
		F(temp, z, distancia, m, M, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K1);

		//stage 2
		produtoEscalarVetor(K1,0.5,tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		F(temp, temp2, distancia, m,M, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K2);

		//stage 3
		produtoEscalarVetor(K2,0.5,tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		F(temp, temp2, distancia, m,M, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K3);

		//stage 4
		produtoEscalarVetor(K3,1.0,tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		F(temp, temp2, distancia, m,M, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,K4);

		produtoEscalarVetor(K1,1.0/6,tamanho,temp);
		somaVetores(z,temp,tamanho,z);
		produtoEscalarVetor(K2,1.0/3,tamanho,temp);
		somaVetores(z,temp,tamanho,z);
		produtoEscalarVetor(K3,1.0/3,tamanho,temp);
		somaVetores(z,temp,tamanho,z);
		produtoEscalarVetor(K4,1.0/6,tamanho,temp);
		somaVetores(z,temp,tamanho,z);

		norma[1]=calculaNorma(z+tamanho/2,tamanho/2);

		if(i%escrever==0)
		{
			sprintf(aux,"%0.15lf ",norma[1]);
			fputs(aux,f[2*corpos+2]);
			fputs("\n",f[2*corpos+2]);
		}

		y[1]=((double)(i-1)/(double)(i))*y[0]+2*log(norma[1]/norma[0]);
		Y[1]=((double)(i-1)*Y[0]+y[1])/(double)i;

		if(i%renormalizacao==0)
		{
			for(j=tamanho/2;j<tamanho;j++)
			{
				z[j]=z[j]/norma[1];
			}
			norma[1]=1;
		}

		if(i%escrever==0)
		{
			escreverArquivos();
		}

		y[0]=y[1];
		Y[0]=Y[1];
		norma[0]=norma[1];
	}

	int minutos=0;
	int segundos=0;
	int horas=0;
	tempo=0;
	int tempo2=0;

	fim=clock();
	total=fim-inicio;
	tempo=(double)total/CLOCKS_PER_SEC;
	tempo2=(int)tempo;
	horas=(tempo2/3600);
	minutos=(tempo2/60)-horas*60;
	segundos=(tempo2 % 60);
	printf("%d h %d min %d s\n",horas, minutos ,segundos);
	printf("%lf \n",tempo);

	sprintf(aux,"%0.10lf ",tempo);
	fputs(aux,f[2*corpos]);
	fputs("\n",f[2*corpos]);

	finalizar();

	return 0;
}

void inicializacao(void)
{
	dimensao=2;
	corpos=3;
	renormalizacao=1;
	T=(int)pow(2.0,n)+1;
	dt=tFinal/(T-1);
	tamanho=corpos*dimensao*2*2;

	int i;

	z =(double *)calloc(tamanho, sizeof(double));
	temp =(double *)calloc(tamanho, sizeof(double));
	temp2 =(double *)calloc(tamanho, sizeof(double));
	dx =(double *)calloc(tamanho, sizeof(double));
	distancia =(double *)calloc(corpos*corpos, sizeof(double));
	Y=(double*)calloc(2,sizeof(double));
	y=(double*)calloc(2,sizeof(double));
	norma=(double*)calloc(2,sizeof(double));
	J=(double *)calloc(pow(tamanho/2,2), sizeof(double));
	K1=(double*)calloc(tamanho,sizeof(double));
	K2=(double*)calloc(tamanho,sizeof(double));
	K3=(double*)calloc(tamanho,sizeof(double));
	K4=(double*)calloc(tamanho,sizeof(double));
	K5=(double*)calloc(tamanho,sizeof(double));
	K6=(double*)calloc(tamanho,sizeof(double));
	aux =(char *)calloc(1000, sizeof(char));

	//delta inicial aleatorio
	srand( (unsigned)time(NULL) );
	for(i=tamanho/2;i<tamanho;i++)
	{
		//delta[i]=rand()/(double)RAND_MAX;
		//delta[i]=delta[i]-0.5;
		z[i]=1.0;
	}

	// norma delta inicial nao normalizado
	norma[0]=calculaNorma(z+tamanho/2,tamanho/2);
	//normalizacao delta inicial
	for(i=tamanho/2;i<tamanho;i++)
	{
		z[i]=z[i]/norma[0];
	}
	imprimir(z+tamanho/2,tamanho/2);
	//norma delta inicial normalizado
	norma[0]=calculaNorma(z+tamanho/2,tamanho/2);

	sprintf(aux,"%0.15lf ",norma[0]);
	fputs(aux,f[2*corpos+2]);
	fputs("\n",f[2*corpos+2]);

	keplerT();

	//velocidade centro de massa ser nula
	for(i=0;i<dimensao;i++)
	{
		z[corpos*dimensao+i]=-z[(corpos+1)*dimensao+i]-z[(corpos+2)*dimensao+i];
	}
}

void keplerT(void)
{
	exemplo="C:\\dados\\";
	mkdir(exemplo);
	exemplo="C:\\dados\\keplerRK4\\";
	mkdir(exemplo);

	M=mSol;
	m=mTerra;
	z[0]=perielioTerra;
	z[3]=mTerra*vPerielioTerra;
}

void criarArquivos(void)
{
	char *raiz,*diretorio,*aux,*barra;

	raiz=(char*)calloc(100,sizeof(char));
	aux=(char*)calloc(100,sizeof(char));
	diretorio=(char*)calloc(100,sizeof(char));

	if(!strcmp(SO,"w"))
	{
		barra="\\";
	}
	if(!strcmp(SO,"l"))
	{
		barra="/";
	}

	strcpy(raiz,exemplo);
	strcat(raiz,barra);
	//mkdir(raiz,S_IRWXU);
	mkdir(raiz);
	sprintf(aux,"%0.0lf",tFinal);
	strcat(raiz,aux);
	sprintf(aux,"%d",n);
	strcat(raiz,barra);
	//mkdir(raiz,S_IRWXU);
	mkdir(raiz);
	strcat(raiz,aux);
	strcat(raiz,barra);
	//mkdir(raiz,S_IRWXU);
	mkdir(raiz);

	if(orbitas)
	{
		strcpy(diretorio,raiz);
		strcat(diretorio,"r.txt");
		f[0] = fopen(diretorio,"w");
		if(f[0] == NULL)
		{
			printf("Erro ao abrir arquivo r!!!\n");
			exit(1);
		}
		strcpy(diretorio,raiz);
		strcat(diretorio,"p.txt");
		f[1] = fopen(diretorio,"w");
		if(f[1] == NULL)
		{
			printf("Erro ao abrir arquivo p!!!\n");
			exit(1);
		}
	}

	strcpy(diretorio,raiz);
	strcat(diretorio,"tempo");
	strcat(diretorio,".txt");
	f[2*corpos] = fopen(diretorio,"w");
	if(f[2*corpos] == NULL)
	{
		printf("Erro ao abrir arquivo5!!!\n");
		exit(1);
	}

	strcpy(diretorio,raiz);
	strcat(diretorio,"megno");
	strcat(diretorio,".txt");
	f[2*corpos+1] = fopen(diretorio,"w");
	if(f[2*corpos+1] == NULL)
	{
		printf("Erro ao abrir arquivo megno!!!\n");
		exit(1);
	}

	strcpy(diretorio,raiz);
	strcat(diretorio,"normaDelta");
	strcat(diretorio,".txt");
	f[2*corpos+2] = fopen(diretorio,"w");
	if(f[2*corpos+2] == NULL)
	{
		printf("Erro ao abrir arquivo normaDelta!!!\n");
		exit(1);
	}

	free(aux);
	free(raiz);
	free(diretorio);
}

void F(double *dx, double *z, double *distancia, double m,double M, int dimensao)
{

	dx[0]=z[2]/m;
	dx[1]=z[3]/m;

	calcularDistancia(z,distancia,dimensao,corpos);

	dx[2]=      -(G*m*M*z[0])/(double)(pow(distancia[0],3));
	dx[3]=      -(G*m*M*z[1])/(double)(pow(distancia[0],3));

	jacobianaK(J,M,m,z[0],z[1],2*dimensao*corpos);

	produtoMatrizVetor(J,z+2*corpos*dimensao,2*corpos*dimensao,dx+corpos*dimensao*2);
}

void calcularDistancia(double *z,double *distancia, int dimensao, int corpos)
{
	int k;
	double temp;

	distancia[0]=0;
	for(k=0;k<dimensao;k++)
	{
		temp=z[k];
		temp=temp*temp;
		distancia[0]=distancia[0]+temp;
	}
	distancia[0]=(double)sqrt((double)distancia[0]);
}

void jacobianaK(double* j, double M,double m,double rx,double ry,int tamanho)
{
	double G2=(double)G;
	double t2 = 1.0/m;
	double t3 = rx*rx;
	double t4 = ry*ry;
	double t5 = t3+t4;
	double t6 = 1.0/(pow(t5,(5.0/2.0)));
	double t7 = G2*M*m*rx*ry*t6*3.0;
	double t8 = 1.0/(pow(t5,(3.0/2.0)));

	j[0*tamanho+2]=t2;
	j[1*tamanho+3]=t2;
	j[2*tamanho+0]=-G2*M*m*t8+G*M*m*t3*t6*3.0;
	j[2*tamanho+1]=t7;
	j[3*tamanho+0]=t7;
	j[3*tamanho+1]=-G2*M*m*t8+G*M*m*t4*t6*3.0;

}

void produtoEscalar(double* v,double* v2,int tamanho, double *resultado)
{
	int i;

	*resultado=0;
	for(i=0;i<tamanho;i++)
	{
		*resultado+=v[i]*v2[i];
	}
}

void somaVetores(double* v,double* v2,int tamanho, double *resultado)
{
	int i;

	for(i=0;i<tamanho;i++)
	{
		resultado[i]=v[i]+v2[i];
	}
}

void produtoEscalarVetor(double* v,double a,int tamanho,double *resultado)
{
	int i;

	for(i=0;i<tamanho;i++)
	{
		resultado[i]=a*v[i];
	}
}

void produtoMatrizVetor(double *J,double* delta,int tamanho, double *resultado)
{
	int i,j;
	double temp;

	for(i=0;i<tamanho;i++)
	{
		temp=J[i*tamanho]*delta[0];
		for(j=1;j<tamanho;j++)
		{
			temp+=J[i*tamanho+j]*delta[j];
		}
		resultado[i]=temp;
	}
}

void escreverArquivos(void)
{
	int j,k;
	char aux[100];

	if(orbitas)
	{
		for(k=0;k<corpos;k++)
		{
			for(j=0;j<dimensao;j++)
			{
				//arquivo txt
				//posicao
				sprintf(aux,"%0.20lf ",z[k*dimensao+j]);
				fputs(aux,f[k]);
				//momento
				sprintf(aux,"%0.20lf ",z[(k+corpos)*dimensao+j]);
				fputs(aux,f[k+corpos]);
			}
			fputs("\n",f[k]);
			fputs("\n",f[k+corpos]);
		}
	}
	sprintf(aux,"%0.15lf ",Y[1]);
	fputs(aux,f[2*corpos+1]);
	fputs("\n",f[2*corpos+1]);
}

void imprimir(double* v,int tamanho)
{
	int i;
	for(i=0;i<tamanho;i++)
	{
		printf("%0.10lf\n",v[i]);
	}
	puts("");
}

double calculaNorma(double *v,int dimensao)
{
	double temp=0;
	int i;

	for(i=0;i<dimensao;i++)
	{
		temp+=v[i]*v[i];
	}
	temp=sqrt(temp);
	return temp;
}

void criarMatrizRK45(void)
{
	int stages=6;
	a=(double*)calloc(stages*stages, sizeof(double));
	b4=(double*)calloc(stages, sizeof(double));
	b5=(double*)calloc(stages, sizeof(double));
	c2=(double*)calloc(stages, sizeof(double));

	a[1*stages+0]=1.0/4;
	a[2*stages+0]=3.0/32;
	a[2*stages+1]=9.0/32;
	a[3*stages+0]=1932.0/2197;
	a[3*stages+1]=-7200.0/2197;
	a[3*stages+2]=7296.0/2197;
	a[4*stages+0]=439.0/216;
	a[4*stages+1]=-8.0;
	a[4*stages+2]=3680.0/513;
	a[4*stages+3]=-845.0/4104;
	a[5*stages+0]=-8.0/27;
	a[5*stages+1]=2.0;
	a[5*stages+2]=-3544.0/2565;
	a[5*stages+3]=1859.0/4104;
	a[5*stages+4]=-11.0/40;

	b4[0]=25.0/216;
	b4[2]=1408.0/2565;
	b4[3]=2197.0/4104;
	b4[4]=-1.0/5;

	b5[0]=16.0/135;
	b5[2]=6656.0/12825;
	b5[3]=28561.0/56430;
	b5[4]=-9.0/50;
	b5[5]=2.0/55;

	//independente
	c2[0]=0.0;
	c2[1]=1.0/4;
	c2[2]=3.0/8;
	c2[3]=12.0/13;
	c2[4]=1.0;
	c2[5]=1.0/2;
}

void finalizar(void)
{
	if(orbitas)
	{
		fclose(f[0]);
		fclose(f[1]);
		fclose(f[2]);
		fclose(f[3]);
		fclose(f[4]);
		fclose(f[5]);
	}
	fclose(f[6]);
	fclose(f[7]);
	fclose(f[8]);

	free(distancia);
	free(a);
	free(b4);
	free(b5);
	free(c2);
	free(K1);
	free(K2);
	free(K3);
	free(K4);
	free(K5);
	free(K6);
	free(J);
	free(temp);
	free(norma);
	free(Y);
	free(y);
	free(z);
}
